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SUMMARY 


Many  modern  compuLiiioruil  fluid  dynomici;  compuier  programs  aft  developed  by  using 
cite  finite  volume  duicreuzdiion  method  It  ha>i  en  excelleni  numencdJ  capability  for  capturing 
changeii  in  conserved  quantities  such  ai  mass,  momenium  and  energy.  In  many  canes,  thermal 
energy  is  irannfened  from  fluids  lo  some  adjacent  solid  mass  Heal  accumulation  in  this  solid 
matter  is  an  imponani  engineering  issue.  To  capture  this  energy  transfer,  it  in  important  to  have 
heat  conduction  algorithms  that  function  well  with  fluid  dynamics  codes.  The  present  work 
tackles  this  problem  by  presenting  an  algorithm  for  solving  the  heat  ei^uation  in  finite  volume 
form.  Although  thin  derivation  is  cast  in  two  dimensions,  it  may  be  readily  generalized  to  three 
dimensions.  Example  problems  are  solved  involving  heal  conduction  within  a  section  of  an 
annular  ring.  Along  the  boundaries  we  enforce  both  Dirichlei  and  Neumann  boundary  conditions. 
For  code  validation,  our  numerical  solutions,  based  upon  the  Ekiuglas-Rachford  ADI  time 
integration  scheme,  are  compared  with  exact  mathematical  solutions. 
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i.»  intr()di;i:ti()n 


Ttte  iranifer  of  heal  been  of  greai  iniere&i  ^viihm  ctte  engineering  and  scientific 
communiiies  for  a  very  long  lime.  Engineering  &>'«tems  art  often  suvepiible  to  ddnuge  from 
extreme  temperaiure.^!,  either  ttoi  or  cold,  so  in  the  course  of  operation,  it  is  imponani  to 
understand  the  temperature  distnbuiion  through  the  material  compnsing  the  system.  Consider  the 
simple  example  of  a  piping  section  that  transfers  fluid  through  a  conventionally  fueled  power 
plant  If  it  leads  from  the  boiler  to  the  turbo-electric  generator,  the  piping  becomes  very  hot 
Since  this  piping  is  connected  to  nearby  thermally  conductive  structures,  we  realize  that  the 
entire  system  has  a  complicated  temperature  distribution.  Temperature  gradients  throughout  the 
system  create  thermal  stress  fields  in  the  piping  material.  Like  its  mechanical  counterparL 
thermal  stress  tends  to  distort  the  size  and  shape  of  a  mechanical  componeni.fi]  Moreover, 
excessive  temperature  gradients  can  cause  serious  component  dar^iage.  If  the  boiler  tubing  is 
suddenly  sub^cted  to  a  torrent  of  cold  water  upon  a  sudden  loss  of  the  heat  source,  the  tubing 
may  crack  or  shatter.  The  piping  may  also  bend  or  change  length  breaking  piping  jomis  These 
deleterious  results  are  ultimately  tied  lo  the  temperature  distribution  in  the  maierial.  Since  society 
relies  upon  ihe  generation  of  energy  and  heal  (either  as  a  direct  or  indirect  product),  the 
importance  of  understanding  heat  transfer  effects  is  made  clear. 

To  provide  clarity  in  the  absence  of  rigor,  we  may  differentiate  between  two  of  the  major 
types  of  heat  transfer  Connective  heat  transfer  involves  the  large  scale  motion  of  material 
carrying  thermal  energy.(2]  The  piping  system  mentioned  above  carries  high  temperature  fluid 
from  a  hot  source  to  a  cooler  beat  sink.  The  motion  of  the  fluid  in  the  pipe  characterizes  this 
transfer  as  being  convective  Conductive  heat  transfer  is  not  the  same;  instead,  it  relies  upon  tiny 
molecular  motions  |J]  (f  you  hold  an  ice  cube  in  your  hands,  heat  conduction  becomes  apparent 
The  ice  is.  for  all  intent,  solid  and  non-moving,  so  no  large  scale  matenal  motion  occurs.  Yet. 
your  hand  becomes  cold  where  it  is  in  contact  with  the  ice,  and  the  ice  begins  to  melt  into  liquid. 
At  the  interface  between  skin  and  ice,  molecules  are  in  contact,  and  the  warmer  skin  molecules 
transfer  kinetic  energy  via  collisions  into  the  colder,  more  lethargic  ice  molecules.  This  process 
1%  governed  by  molecular  diffusion,  not  by  large  scale  material  motion.  Tlie  work  that  follows 
serves  to  model  heat  conduction  in  a  careful  way  with  a  flexible  set  of  numerical  algorithms 
suitable  for  use  in  solving  problems  involving  non-Canesian  geometnes. 

Standard  numerical  techniques  for  modeling  heal  conduction  are  really  developed  for 
Cartesian  coordinates,  but  the  vast  ma^rity  of  engineering  problems  have  non-Cartesian 
geometries.  A  circular  walled  pipe  is  a  good  example  We  can  only  apply  Cartesian  algorithms  to 
this  problem  if  we  apply  special  computer  coding  lecluiques  such  as  ti>  using  a  system  of 
interpolations  at  die  boundary  to  remap  die  cylindrical  problem  geometry  onto  a  background 
Cartesian  gnd,  or  (in  representing  the  body  as  complicated  interface  on  the  Cartesian  grid.  This 
interface  is  effectively  tracked  on  the  Cartesian  gnd.  and  the  physics  of  the  problem  can  dien  be 
stored  in  reference  to  Cartesian  grid  Approach  {i)  is  quite  old,  and  it  is  rife  with  error  due  to  the 
need  for  an  interpolation  remapping  scheme.  |4]  This  scheme,  in  general,  ha.s  no  relevance  to  the 
physics  of  die  problem.  Near  oddly  shaped  boundaries,  the  use  of  interpolation  also  requires  a 
great  deal  of  bookkeeping  to  track  the  indicial  limns  of  the  boundary  on  die  grid.  Approach  {li)  is 
quite  valuable,  yet  it  is  still  a  matter  of  contemporary  research  in  numerical  algorithms  |5] 
Modern  interface  tracking  methods  are  very  complicated,  but  they  can  apply  to  a  wide  class  of 
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problems.  In  fact,  you  can  treai  mulii-physics  problems  by  using  this  approach.  Ho^vever,  ii 
dimply  noi  necessary  to  reson  to  ihi.«  much  nun^erical  machinery  in  order  lo  develop  a  rodu^i  and 
accurate  compuuiiotul  solver  for  heal  conduction.  Instead,  the  finite  volume  discretization 
technique  may  be  applied  in  conjunction  with  a  generalized  coordinate  transformation  in  order  to 
map  the  governing  equation  onto  a  system  of  body-fitted  coordinates.  The  transformed  system  is 
readily  solvable  by  conventional  discretization  techniques.  In  fact,  it  is  desirable  that  our  time 
and  space  discretization  equations  look  as  much  like  dieir  clasiical  Cartesian  counterparts  as  is 
possible.  This  widely  applicable  approach  lies  at  the  foundation  of  modem  numerical  disciplines 
like  computational  fluid  dynamics  ICFD)  and  computational  electromagnetics  <CEM>  Like  these 
disciplines,  computational  heat  transfer  relies  upon  the  continuum  mechanics  interpretation  of 
heal  flow  through  governing  partial  differential  equations. 

1.1  Efqualion  for  Heal  Conduction 

Lei  K  be  the  temperature  ai  a  point  in  space  within  a  mass  of  material.  Generally, 
temperature  is  measured  on  an  absolute  scale,  but  relative  scales  are  also  acceptable  here.  The 
equation  representing  the  time  dependent  conduction  of  heat  (without  the  presence  of  internal 
heal  sources),  as  represented  by  the  temperature,  is  written  as  follows  |6] 

=  0  (1) 

ln(l>,  UBU(j;y>s,r),  where  x,  s'  and  ;  are  the  Cartesian  coordinates  of  the  point:  the  time  is 
denoted  as  / .  The  parameter  a  the  coefficient  of  thermal  diffusiviiy  defined  as 


(2) 


where  k  is  the  coefficient  of  thermal  conductivity  for  the  material  (taken  as  constant  in  this 
case),  while  p  and  C«  are  the  density  and  specific  heat  capacity  of  the  material  For  problems 
involving  a  spatially  dependent  thermal  conductivity,  ( I  >  has  a  slightly  differeni  form,  but  our 
numerical  algorithms  function  just  as  well  for  that  case.  Equation  (1)  has  a  number  of  names 
such  as  ihe  heal  equation,  the  ditfusion  equation  and  Fourier's  equation  More  imponanily,  it  is 
classified  as  a  linear,  parabolic  paniaJ  differential  equation.  This  classification  is  used  to  detail  a 
set  of  maihemaiical  techniques  that  may  be  used  to  obtain  "exact”  solutions  for  ( I ).  Later  in  this 
report,  we  apply  some  of  these  techniques  to  build  validation  solutions  for  our  numerical 
algorithms.  Equation  (1)  also  contains  the  Laplace  operator  defined  as 


dx  dy  dt 


(3) 


m  three  Cartesian  dimensions.  For  problems  ca.st  in  fewer  dimensions,  we  simply  delete  from  (3) 

the  derivative  for  the  missing  coordinaie(s>  The  Laplace  operator  can  be  expressed  in 

generalized  (non-Canesian)  coordinates,  bui  its  mathematical  form  is  greatly  complicated  by  this 

process  The  heat  equation  (or  diffusion  operator)  tends  to  add  a  great  deal  of  smoothness  lo  an 
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iniiial  lemperuure  distribution  (or  to  (he  distribution  of  some  other  scoldi  quantity).  Aldiough  (3) 
1%  relatively  easy  lo  dii^cretize  by  using  fuuie  differences,  its  form  in  generaliaed  coordinaie^  is 
not.  Later,  we  discuss  the  application  of  (he  finite  volume  discretization  method,  a  technique  we 
can  use  to  overcome  this  difficulty. 

The  heat  equation  is  normally  solved  as  an  initial-boundary  value  problem  (IBVF) 
meaning  (hat  (he  equation  must  be  accompanied  by  (i)  a  domain,  (ii)  a  set  of  conditions 
specifying  either  the  solution  or  ns  derivatives  along  the  boundary  and  (iii>  initial  conditions 
specifying  (he  solution  at  every  point  in  the  domain  at  the  starting  time  for  (he  problem.  In  one 
dimension,  an  example  of  a  classical  IBVP  for  the  heat  equation  is  shown  below. 

—  0<x<l.  i>0 

d> 

M0./)s0  uCl,/)sO  (4) 

uUQUfix) 

In  (4),  the  diermal  diffusiviiy  is  set  equal  to  one.  The  time  and  space  domains  are  clearly 
described  in  (4>,  and  the  boundary  conditions  are  provided  on  the  second  line.  Note  that  the 
boundary  consists  of  exactly  two  space  points  for  the  problem.  The  initial  conditions  are 
delineated  on  the  last  line.  For  a  piecewise  continuous  function  /{x),  (his  problem  is  easily 
solved  by  separating  variables  and  using  Fourier  series.  The  exact  solution  is 


o; 

a{x,  0  s  expi-^j  V*/ )  8in(«rr  x)  (5) 

where 

1 

4.  «2j/U)sin(/iffx)r£r  (6) 

0 

Exact  solutions  may  be  obtained  for  the  heat  equation  in  many  cases,  but  only  for  problems  cast 
in  simple  orthogonal  coordinate  systems  <e.g.,  Cartesian.  Cylindrical,  Spherical).  For  other 
coordinate  systems,  we  must  resort  to  numerical  solution  techniques.  This  topic  is  the  sublet  of 
(he  remainder  of  this  report. 

1.2  A  Map  for  (he  Remainder  of  IliU  Report 

In  Section  2,  a  discussion  of  the  finite  volume  method  is  presented  along  with  complete 
derivations  required  for  discretizing  the  heat  equation  on  non-Cartesian  domains  in  two 
dimensions.  Although  this  exposition  is  limited  to  two  dimensions,  the  method  is  easily  extended 
to  three  dimensions  and  for  problems  involving  non-constani  thermal  conductivities.  To  assist  in 
discretizing  curvilinear  bodies,  a  local  coordinate  transformation  is  introduced  and  used  in 
transforming  derivatives  With  the  spatial  discretization  scheme  intact,  we  express  finite  volume 
interface  quantities  in  a  form  suitable  for  use  m  alternating  direction  implicit  (ADI)  and  explicit 
time-integration  schemes  In  fact,  a  finite  volume  interpreiaiion  of  the  Douglas-Rachford  sc^me 
IS  presented.  Section  3  contains  a  series  of  test  cases  for  the  numerical  algorithm.  The  basic 
numerical  concept  is  demonstrated  first  for  rectangular  coordinates  in  solving  the  asymmetrical 
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heuing  pfoblem.  Although  this  test  ca&e  is  i^imple,  it  posse^es  resl  uuliry  m  testing  the 
numenca)  scheme  Secondly,  we  solve  the  heat  ei^uaiion  on  a  90’  sector  of  an  annular  ring 
undergoing  asymmetrical  heating  m  the  radial  direction.  This  problem  is  characterized  by 
azimuthal  (or  angubr)  symmetry,  yet  it  provides  a  thorough  teiit  of  the  time-marching  i^cheme. 
For  the  final  leiit  case,  ^ve  again  employ  the  annular  nng.  but  without  azimuthal  symmetry. 
Asymmetric  heating  is  applieil  in  both  the  radial  and  angular  directions,  particularly  along  the 
outer  boundary  of  the  ring  sector.  For  each  test  case,  the  numerical  solution  is  compared  with  an 
exact  mathemaiioal  solution  determined  by  Fourier  series  analysis.  Good  agreement  with  the 
exact  solution  is  demonstrated  for  each  test  case,  validating  our  numerical  methods  Section  4 
contains  concluding  remarks  that  pertain  largely  to  the  efficacy  of  this  numerical  approach  to 
heal  conduction.  Cautionary  remarks  are  made  on  the  treatment  of  non-orthogonal  curvilinear 
coordinate  systems.  Large  numerical  discretizaijon  errors  can  occur  in  these  coordinate  systems 
if  the  discretization  is  not  conducted  with  exceptional  care.  Time  dependent  boundary  conditions 
are  also  discussed  since  the  temperature  (or  heat  flux)  is  forced  to  change  in  time  along  the 
boundary  of  many  engineering  devices.  Near  the  end  of  this  report,  appendices  are  included  to 
describe  the  derivation  of  the  exact  solutions  These  solutions  are  thoroughly  discussed  because 
of  their  importance  in  the  validation  process. 
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2.0  METHODS,  ASSUMPTIONS  AND  PROCEDURES 


Ttte  etposiiion  ihai  follows  conidin^  ihe  nidittemaiiciil  machinei^  required  u>  trdn>;form 
the  heal  equation  from  differential  form  into  a  set  of  algebraic  equations  that  may  easily  be  coded 
in  a  computer  language  Diligence  is  applied  in  maintaining  the  derivation  as  rigorous  is  ai 
possible  without  xacrificing  practicality.  At  the  end  of  the  day,  a  report  that  lacks  implementation 
details  is  of  little  use  to  the  practicing  engineer  or  scientist.  Naturally,  we  avoid  mathematical 
proofs  involving  concepts  such  as  convergence  and  stability,  and  this  approach  is  reasonable 
since  we  use  legacy  numerical  schemes  (admittedly  with  some  required  and  well  documented 
changes)  whenever  possible  Details  associated  with  the  mechanics  of  our  finite  volume 
discretization  we  endeavor  to  explain  in  detail 

2.1  Finite  Volume  Discretisation 


Many  scientists  are  familiar  with  the  finite  difference  method  for  discretizing  differential 
equations.  In  fact,  the  Newton  divided  difference  is  the  numerical  analog  for  a  denvative  |7] 
Interestingly  enough,  the  finite  volume  method  relies  upon  integration  via  the  Fundamental 
Theorem  of  Calculus  for  the  simple  one-dimensional  example  shown  below.  Suppose  we  wish  to 
evaluate  dyidx  at  some  point  x,  Let  us  integrate  this  derivative  on  the  interval  . 


Also,  we  can  envidon  the  concept  of  an  averaged  value  for  ttytdx  defined  at  x, .  In  this  sense, 
we  have  that 


•,j.»  I*, 


Hence. 


V2 

dx  it 


<9) 


where  The  expression  (8)  for  the  average  derivative  matches  the  Newton 

divided  difference  formula,  so  for  uniform  one-dimensional  meshes,  the  finite  volume  and  finite 
difference  discretization  procedures  render  the  same  results.  It  is  insinictive  also  to  note  that  the 
finite  volume  method  introduces  the  concept  of  "averaged”  values.  More  imponaniiy,  the  tiniie 
volume  procedure  has  even  greater  utility  in  higher  dimensions 


Consider  an  example  cast  in  a  plane  as  is  shown  in  Figure  I.  For  simplicity,  we  have 
chosen  noi  to  show  a  complete  body  (or  domain)  geometry.  Instead,  we  have  outlined  a  single 
finite  volume  "cell"  m  dashed  lines.  Not  that  ihe  sides  of  ihe  cell  are  not  aligned  with  any  set  of 
Cartesian  coordinates.  The  cell  is  identified  by  indicial  coordinates  {Kj).  In  the  discussions  that 
follow,  the  vertex  centered  finite  volume  method  is  applied.  That  is  to  say,  die  center  of  cell 
coincides  with  a  grid  point  possessing  the  same  indices.  The  cell  centered  finite  volume  is 
equally  valid  and  invokes  the  so-called  "northeast''  cell  numbering  convention.  According  to  this 
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convention,  ihe  grid  pomi  loisigns  Us  indices  lu  the  coll  center  to  its  noithcasL  The  sA^igmng  grid 
poiAi  is  circled  in  Figure  I.  These  indicial  coofdituie«  are  associated  with  o  system  of  local 

\  > 


^  s 

Figure  I '  ConTiguraiion  of  a  VertOK  Centered  QuadnlateraJ  2D  Finite  Volume  Cell 


coordinates  1^.  t))  sltowe  m  the  figure  The  local  coordinates  cuo  be  identified  with  direcdonal 
* 

vectors  |  atid  ;; .  The  directions  of  these  isso  vectors  cluuige  with  each  point  in  the  grid  unlike  a 
Cartesian  system. 


iH 


where 


cx  5v 


ac*  ay* 

7  being  the  Grudieni,  a  commonly  used  vector  difTercntiaJ  operator. 


(10) 

(11) 


To  begin  die  finite  volume  discretization  of  the  2D  heat  et|uation  above,  we  invoke  the 
Divergence  Theorem,  a  higher  dlmenniunal  version  of  the  Fundamental  Theorem  of  Culculus.[8) 
Integrate  (10)  over  the  interior  volume  of  the  cell  Note  that  in  2D,  volume  L<  interpreted  as  area, 
and  it  follows  that  the  2D  analog  of  surface  area  (for  die  cell  sides)  is  length.  Accordingly,  the 
quudriluientl  cell  shown  in  Figure  1  has  four  sides.  labeled  I  through  4  Each  cell  side  m  (or 
Interface)  has  length  m  =  I...  ,4  We  have  that 

J^Ji'  =  Jv  Vutfv  =  ^  Vd'iidr  (12) 

The  nglit  most  term  in  (12)  is  has  been  produced  by  applicalion  of  the  Kvergence  Theorem.  Tbe 
overall  volume  (area)  integral  bait  bees  iraosforreed  into  as  integral  over  the  closed  surface  (side 
leogihi)  of  the  cell.  The  symbol  n  represents  die  outer  pointing  unit  normal  vector  defined  at 
eoch  point  on  the  cell  interface.  Since  the  cell  has  four  straight  sides,  r)  takes  on  four  distinct 
vector  vulue.s.  The  left  hand  side  of  (12)  can  also  be  rewiiiien  For  our  purpose!;,  the  grid  does 
not  move,  so  the  order  of  differentiaiion  and  iotegretion  can  be  interchanged,  i.e . 

=  L\ad\>  (131 

J  A  di  * 


Oemei/nON  a.  A9<6»M  Io<  pucnc  r«nt».asl;»uMf«ur^dM  9e''AB'IV'PA  Aesreua  M  CMva  AaeABW  SOKMUOd. 
OaMSAcriSmO  9 


Recall  ihai  UBu(xy.r),  &o  ihe  process  of  iniegraiing  u  over  cell  volume  defmes  a  tie  fin,  iv 
average  C  in  thai 

=  uV  <I4) 

By  «ubsiiiuiing<IJ>atKl  (14)  into  (12).  we  have  (hai 


Note  Out  because  of  Ok  process  of  spatial  integration,  Q, ,  does  not  depetvj  on  any  particular 


point  inside  of  the  cell:  in  the  cell,  it  ui  solely  a  function  of  time  Hence,  a  total  time  derivative 
has  replaced  the  partial  denvauve.  The  symbol  V  ,  represents  the  total  volume  (or  area,  in  this 

case)  for  die  cell  Since  the  cell  boundary  is  represented  by  a  sei^uence  of  line  segments,  we  can 
simplify  the  right  hand  side  of  <  15).  (n  fact,  ihe  integral  over  the  dosed  boundary  can  be  replaced 
by  a  summation  of  integrals  along  directed  line  segments,  ik.. 


At  this  point,  the  finite  volume  method  applKS  an  approximation.  If  the  function  u  was  known  as 
a  continuum  in  an  open  neighbor  of  the  cell  sides,  the  integration  in  ( 16)  could  be  performed  in 
an  exact  manner.  UnfortunaKly,  u  is  usually  known  only  at  a  discrete  number  of  points. 
Moreover,  in  order  for  (16)  to  form  a  proper  numerical  method,  we  must  moruiruu  from 
the  locus  of  S .  This  approximation  is  made  at  the  cenKr  of  each  side  of  the  cell.  That  is  to  say, 
Vu  is  estimated  at  the  center  point  of  each  cell  side  then  multiplied  by  the  side  length.  Likewise, 
the  normal  vector  ri  is  calculated  at  the  same  point.  Equation  (16)  then  becomes 


du 


di 


ViT-n.L 


(17) 


In  (17),  geometrical  quantities  such  as  cell  volume,  side  normal  vectors  and  lengths  are  easily 
calculated  based  upon  the  grid  point  locations  while  quantities  based  upon  S  require  a  more 
considered  approach  When  a  diicreti2ed  form  for  Vi7  is  used  (a  sublet  addressed  later),  one 
derivative  remains  in  this  expression.  For  that  reason,  ( I7>  is  referred  to  as  a  semi-discreie  form 
This  form  is  very  useful  since  it  is  amenable  to  a  number  of  different  temporal  discretuation 
schemes.  Still,  we  must  address  differentiation  on  the  curvilinear  grid  before  we  can  completely 
formulate  a  discretized  version  of  . 


2.2  Traiirtfortnalion  Mathematics  and  Local  Coordinates 

The  use  of  Cartesian  discretization  is  simply  impraciical  for  most  real  world  engineering 
problems.  Even  when  appiKd  along  with  a  modem  interface  tracking  method  for  oddly  shaped 
bodies,  interesting  questions  arise  regarding  the  accuracy  of  derivatives  computed  near  the 
boundary.  For  instance,  if  a  thermal  insulator  is  installed  on  the  curved  boundary  of  a  maKnal, 
physics  dictates  that  we  enforce  boundary  conditions  on  heat  flux  in  the  direction  of  the  normal 
vector  ai  every  point  on  tlie  boundary.  A  Cartesian  approximation  of  this  process  enforces  this 
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buimdiiry  condition  by  re&olving  il)e  heoi  fluK  vector  along  two  Caite^ian  coonluuies,  T1u» 
approach  workc,  boi  the  absolute  reliance  upon  fictitious  coordinate  systems  and  polygonal  cut 
cells  tends  lo  leave  the  careful  scieoliii  with  a  feeling  of  unease.  Instead,  let  us  consider  onoiber 
approach  to  the  problem. 

Through  ibe  use  of  calculus,  we  cun  construct  a  transforreatioo  that  maps  the  body  and  its 
boundary  onto  a  Cunesiun-like  computiitioniil  space.  In  thi«  space,  it  is  easy  to  compute 
differences  along  smooth  coordinate  lines  that  conform  to  the  shape  of  the  boundary.  If  the 
curvibneor  coordinates  ore  nearly  orthogonal,  then  derlvtiiives  calculated  with  the  use  of  these 
coordinates  take  on  simple,  easily  discretkoble  forms.  To  support  the  2D  example  shown  in 
Figure  2,  consider  a  iransfumuijon  T  from  the  1^.  r?)  curvilmeor  coordinate  system  to  the 
(i’.y)  Caiiesiao  system,  i£., 

1181 


Rgure  2'  A  Body  Cooformmg  Coordinate  System  fur  a  Seciioo  of  on  Annular  Ring 


The  Jacobian  matrix  for  this  transformation  is  wrirten  os 


The  inverse  function  theorem  guaruntees  that  if  IJkO.  then  the  inverse  trunsrornution 
exists  such  that  7"' ;  (x.  y}^  .17)  .This  transformation  has  Jacobian  moirix  J*'  given  by 


(191 

7..1 


V. 

i,' 

1 

J  V 

»7, 

"iTi 

where 


(20J 


ljl=x.y.-x^y^  (21) 

Note  that  each  point  assigned  within  the  gnd  {a  simple  polar  grid.  In  this  case)  has  a 

corresponding  set  of  (x:.  y)  coordinaieA.  Moreover,  the  body  conforming  grid  has  real 
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^vanL)ge&.  Noie  iJui  the  ^  azimuthal  coofduuie  curve  n  paran^eienzed  by  the  index  i  while 
itte  radidJ  coordiiuie  1/  is  indexed  by  y.  As  a  result  ihe  compuuiioruil  ^id  poini  (^, .r;,) 
direcily  coincideji  with  the  point  physical  space  Using  ihii;  framework, 

derivativeii  such  as  Xg  and  are  easily  calculated  at  grid  point  (^,>7,)  (or  U,j)  in  sltorthand 
notation),  i.e.. 


^  4t+1.y)-J(t-1.y)  . 


V 


y(Ay+l)-y(/>;-i) 

2 


(22) 


As  you  may  note,  the  divided  differences  in  <22>  contain  no  obvious  measurement  of  die  change 
in  ^  or  17 .  This  fact  is  a  testament  to  the  power  of  this  method.  The  mesh  for  the  computational 
plane  is  uniform  dy  definition.  Therefore,  since  both  ^  and  rj  are  elements  of  the  natural 
numbers.  ^=Ar7=1.  Since  the  denominators  in  (22)  are  respectively,  2A,^  and  2A17.  the 
estimates  for  these  derivatives  logically  follow.  Formulas  such  as  those  shown  in  (22),  can  then 
be  used  in  (20)  to  produce  derivatives  such  as  and  rj.  a&soaated  with  7^' .  and  the  set  of 


transformations  has  been  completely  described. 


Our  segue  into  coordinate  transformations  may  have  seemed  to  wander  off  topic,  so  it  li 
now  time  to  demonitrate  its  utility  for  the  heat  conduction  problem  Recall  that  in  ( I7>  the  finite 
volume  scheme  requires  that  be  evaluated.  The  gradient  operator  (II)  is  difficult  to 

evaluate  for  the  boundaries  shown  in  Figure  2,  so  let  us  transform  this  operator  to  the  body  fitted 
(^.17)  coordinate  system.  By  using  die  chain  rule,  we  obtain 

(23) 

[dxd4  dxdrjl  \dyd4  dydrj) 

or  in  our  compact  notation. 


(24) 

Although  we  omitted  the  averaging  overbar  in  (23)  and  (24),  the  gradient  is  evaluated  in  terms  of 
cell  averaged  values  for  u.  Moreover,  for  the  finite  volume  method,  each  cell  is  assumed  to 
possess  a  single,  ^'averaged'  value  of  ir  We  assume  that  there  is  no  internal  variation  for  u 
inside  of  the  cell  Other  numerical  methods  do  permit  such  variation,  (e.g.,  the  Discontinuous 
Galerfcin  Method  [9]),  but  those  metliods  are  beyond  the  scope  of  this  report.  Furthermore,  we 
just  don't  need  that  much  horsepower  Referring  to  Figure  1,  it  is  important  to  realize  that  the 
gradient  is  to  be  evaluated  at  the  center  of  each  cell  side  (numbered  I  through  4).  Also,  according 
to  (17),  the  normal  vector  for  the  cell  side  also  participates.  Since  we  have  conveniently 
employed  a  cell  configuration  with  rectilinear  sides,  die  normal  vector  remains  die  same  along 
the  chosen  side.  Hence,  we  can  evaluate  ri  .  an  expression  for  heat  flux,  as  follows 


Vtf.  fl  =  4.  (^.  +>7.  ir,)+^,  +7.  i(^) 

where  the  unit  normal  vector  n  is  written  as 


(25) 
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(26) 


ns,n,*)n, 

A  icore  compuialjonally  prdciical  form  li 


^u.h  =  (A,‘,+A,‘,)u.+{h,  rj.+n,  rjju^ 

(27) 

^u>n=$.  u,+Au_ 

(28) 

Vp’bere 


The  forniN  shown  in  (28)  uiJ  (29)  ore  qiuie  useful  when  developing  lime-nutrching  kchenies  for 
the  hmi  M|iBtijon.  Aguin.  these  equations  ere  evuludieJ  ai  the  center  point  for  each  cell  side  to 
assist  in  the  heat  Bux  discretiuiiun 


23  Spatial  DUcrellzatioD  of  Heat  Flux 


A  purpose  of  the  senu-discreie  funn  is  to  separate  the  procedures  for  spatial  and  temporal 
disctciiration  allowing  u  choice  of  disctciixniion  methods.  In  this  section,  we  present  the 
spatially  discreti2ed  forms  for  (28)  and  (29).  along  with  the  cell  side  length,  core  features  of  (17). 
As  a  visual  aid.  Figure  3  is  employed  along  with  Figure  2, 


Ir) 


-Ut 


•Cl 


9 


Figure  3;  Grid  Array  for  Vertex  Centered  Finite  Volume  Cell  (ij) 


(n  Figure  3,  the  comers  of  the  cell  ore  indicated  by  red  X's.  The  cell  sides  are  drawn  in  as  dashed 
liues,  ood  gnd  node  (ij)  is  situated  at  the  center  of  cell  (t^).  The  surroundiiig  grid  points  needed 
for  various  discreLzairon  procedures  are  also  included  in  this  diagram  Cell  sides  1  through  4  ore 
labeled  in  blue  For  purposes  of  bookkeeping,  in  the  computational  space,  the  cell  corner  points 
have  hiilf-i30de  indices.  That  is  to  say.  the  upper  left  hand  comer  has  gnd  indices 
the  remaining  comers  are  indexed  accordiogly  Similarly,  the  center  point  for  side  1  con  be 
thought  of  as  having  gnd  indices  .  For  the  vertex  centered  finite  volume  mesh,  the 

cenier  of  cell  (ij)  has  die  coordinaies  (x,  ^.v, ,)  or  (9, >77,).  The  heat  flux  term  requires 
discreiizaiion  on  each  side  of  the  fuuie  volume  cell 
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TtK  ttedi  flux  dL«creU2aiion  for  side  I  is  conducted  iis  follows.  Ai  the  ceaier  point  of  ihe 
«ide,  we  have  ihai 


TtK  form  of  ^ii, ,  depends  on  whether  we  are  employing  an  explicit  or  an  implicit  temporal 
div;reiizaiion  vheme  In  either  case,  the  form  of  (5„u  is  simple,  but  it  is  to  be  diseased  later  in 
this  report  For  now,  we  are  interested  in  the  calculation  of  and  .  These  equations  are  given 
by  (29),  but  the  transformation  derivatives  require  the  following  formulas 

<3I> 

=  (32) 

=i  ('',.,^1  f  )  (33) 

'fl = T  (y.1.^1  -  + y,.«i  -  y,  ,.i )  (34) 

The  side  length  and  normal  vector  components  may  be  calculated  based  upon  the  locations  for 
the  corner  points  for  side  I  The  tangent  vector  directed  along  side  1  is 


^  ~  ^yi2.,.<i2  — '  +  y  (35) 

where  the  notation  i  is  simply  shorthand  representing  the  (xy>  coordinates  of  the  comer  point 
under  consideration  More  specifically. 


®  y(_V*,;.V8~y(_V*,;_V8 


(36) 

(37) 


Care  must  be  exercised  when  calculating  the  normal  vector.  The  Cartesian  coordinate  system 
(x,y.:)  IS  right  handed,  so  it  important  to  ensure  that  the  is  also  right  handed.  Note  that 

we  have  chosen  :  as  our  third  generalized  coordinate  since  the  problem  is  confined  to  ilie  plane. 
Having  observed  this  caution,  the  normal  vector  is  given  by 


(38) 

where  A  is  the  unit  vector  in  the  ;  direction.  Hence,  the  unit  normal  vector  components  are 


where  the  cell  side  length  I  is  calculated  by  using  the  equation 


(39) 


(40) 

Equations  (31)  through  (34)  are  used  in  (20)  and  (21)  permitting  the  computation  of  ,  rj, 
and  r}, .  Final  ly,  by  using  (29),  and  can  be  easily  calculated  for  cel  I  side  1 . 


OfiTflBJTtON  A.  A9<»/«d  IO<  puUiC  r«le«t».aitirbuMA  untmnM  Se''ABW.^  AS&ABW  2010-Oa(M. 


Now  con&ider  side  2  as  is  shown  ia  Figure  3.  T1)e  heal  fluv  lerm  in  (17)  muy  be  wrinen  ds 


The  cakuUuon  of  6^  and  is  now  centered  on  the  midpoini  of  «ide  2.  Tran>;fonnation 
denvative.«  and  other  parameters  can  be  evaluated  with  the  help  of  the  formulas  below. 


•'V  =T(f  -•'-1.  > 

(42) 

y,  =  -X  1.,) 

(43) 

(44) 

(45) 

The  langeni  vector  i  s 

A  A 

^  *  •*— 1/S  j.**  ^  r,  +  /r. 

(46) 

The  nom)al  vector  may  be  calculated  as 

nsitxr 

(47) 

atvj  the  umi  vector  is  given  by 

f.  -  f, 
n  sx — ;  As  — 

‘  /  •  / 

(48) 

where,  as  before,  /sir’ll.  The  iransforntaiion  parameters  leading  to  6, 
calculated  by  using  the  equations  given  in  Section  2.2. 

and  may  be 

For  side  3.  the  expression  for  the  heat  flux  may  be  written  as 

(49) 

The  midpoint  of  side  3  is  the  center  for  the  calculation  of  $  and  the  transformation 

parameters.  Tlie  following  formulas  are  helpful  in  die  calculations. 

(50) 

(51) 

(52) 

= T  (>',.^1  -  y.  ,.i + -  y^y^x ) 

(53) 

In  our  notaiional  convention,  the  tangent  vector  for  side  3  is  calculated  as 

A 

(54) 

A  vector  croei  product  provides  the  ai^sociaied  t>orindl  vector,  x.e.. 

• 

nsrxit 

(55) 
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Hence. 


(56) 


where  I  ^1  f  II.  The  remdjnmg  paran^eiers  occurring  in  (49)  are  calculaied  in  ihe  same  manner  a£ 
that  applied  for  side&  1  and  2. 


The  heal  flux  for  side  4  may  be  calculated  m  ihe  same  way.  Narurally.  caution  must  be 
taken  when  aligning  the  oneniaiion  of  vectors.  The  formula  for  ihe  heat  flux  is 


(57) 


and  the  traniformaiion  parameters  may  be  with  the  aisistance  of  the  following  formula.«. 


.  =  7  (•'..i,  -  -  •'.-t.-i)  <58) 

>: = T  O'-v,  -  y,.i,  ^  -  y^i.,.1 )  (59) 

=•*-.; -•'-M 

,'.  =  y.,-ysH 

TM  vector  for  ^ide  4  may  be  calculated  ai  follows. 


The  disociaied  normal  vector  is  given  by  the  cross  product  nsfyk:  normaliaed,  its 
components  are 


The  above  relations  permii  the  calculation  of  ihe  heat  flux  passing  through  each  of  the  four  sides 
of  the  quadrilateral  finite  volume  cell  illustrated  in  Figure  3.  From  the  standpoint  of  geometry, 
the  formulation  for  heat  flux  is  complete.  SuH.  we  have  not  completely  described  the 
contribution  of  average  cell  temperaiure  IT .  through  terms  ,  to  ihe  heat  flux  formulas  <J0>. 
(41).  (49)  and  (S7).  The  structure  of  this  term  is  dictated  by  the  time  integration  scheme 

2.4  E^xpllclt  niiM  liiteeratlnn 

One  should  not  make  light  of  semi-discreie  formulas  such  ai  (17).  This  form  effectively 
decouples  the  temporal  and  spatial  discreiizaiion  schemes  allowing  a  greai  deal  of  flexibiliiy  in 
algorithm  selection.  The  time  integration  scheme  advances  the  numerical  solution  for  one  instant 
in  lime  to  the  next.  When  using  an  explicii  time  iniegraiion  scheme,  the  numerical  solution  for  u 
«i  the  new  time  level  depends  only  on  the  values  of  i7  at  ihe  previous  time  level.  Explicit 
schemes  can  be  very  accurate.  Put  they  can  also  be  very  difficuli  to  siabiliae.  That  is  to  say,  the 
iruncaiion  error  may  grow  out  of  control  if  the  lime  step  is  loo  large.  For  certain  prodlems, 
explicit  schemes  can  require  prohidiuvely  small  ome  steps  |4]  Implicit  time  integration 
algorithms  are  structurally  different.  For  this  type  of  procedure,  the  numerical  solution  depends 
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upon  (he  values  of  IT  gi  both  ihe  old  and  new  ume  siep«  An  implicit  scheme  generally  requires 
solving  a  system  of  linear  equations.  This  process  can  become  quiie  complex  and  lends  lo  inject 
numerical  viscosity  inio  the  iniegraiion  procedure  On  d)e  other  hand,  explicit  schemes,  even 
those  composed  of  muluple  steps,  are  far  easier  to  program  and  evaluaie.  However,  they  possess 
little  numerical  viscosity  and  are  particularly  unforgiving  of  programming  errors  For  these 
reasons,  explicit  schemes  are  ideal  for  debugging  the  “physics”  represented  by  the  discretized 
equations. 

Our  explicit  time  integration  scheme  as  based  upon  the  modified  Euler  method.|IO] 
Equation  (17)  can  be  written  in  simple  notation,  i  e., 


Js 


ill 


/(ff 


,) 


where  m  is  an  integer  that  is  selected  “cover"  the  extent  of  the  difference  stencil  for  B,  ^ .  This 
equation  can  be  placed  in  modified  Euler  form,  a  Runge-Kuita  scheme  consisting  of  two  steps,  as 
follows 


erf, 

where  >i  is  the  time  step  counter,  i.e,  17'~'  represents  the  solution  at  the  new  time  level.  This 
explicit  scheme  is  formally  of  second  order  in  time,  and  especially  for  testing  spatial  integration 
schemes,  functions  quite  well.  Of  course,  if  higher  order  temporal  accuracy  is  needed,  the 
classical  fourth  order  Runge-Kutta  scheme  can  be  employed.  1 10]  Approximately  twice  as  many 
floating  point  operations  are  required  for  implementing  the  fourth  order  scheme  than  are  required 
for  implementing  (65). 

Numerically,  the  right  hand  side  of  (M)  must  be  evaluated,  in  part,  through  the  use  of  the 
heal  flux  equations  (30),  (41),  (49)  and  (57).  In  these  equations,  (he  terms  and  <(^17  require 
description.  Equation  (30)  represents  the  heat  flux  through  cell  side  I  located  at  computational 
space  coordinates  (/+1/2, )) .  This  interface  is  oriented  in  the  direction  of  ^  (or  J)  coordinate 
lines,  so  we  cannot  build  an  i/  directional  divided  difference  at  this  location.  Instead,  it  is 
necessary  to  average  difference  expressions  evaluated  at  ( /  - 1  y )  and  {i.J  }  Thus. 

“I,  .  V*  ,  *  3  (S,  .1.^1  -  ,.1  -5,  j  <66) 

For  cell  side  3  located  at  computational  space  coordinates  ( /  -f  1  /  2,  jX  is  evaluated  for  (4 1  > 
in  a  similar  manner,  i  e.. 


V  i-v*,  “  > 


(67) 
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To  cwnpuie  the  heal  flux  located  ai  (r,  y -1-1/2  >  for  &ide  2,  (49)  ro^mres  a  detoileil  expre£&ion  foe 
.  Here,  itte  direct  compuLiiion  of  a  divided  difference  along  die  ^  coordinate  i&  not  possible 
For  ihe  explicit  scheme,  averaging  can  be  used  to  obtain  this  term,  i.e.. 


,.irs  *  *  ^ 

Along  side  4,  located  at  (/,  y>1/2>.  we  employ  a  jiimilar  exprei^sion  for  the  divided  difference 
along  the  cell  side 

*  3  +  "-,^1  -  1  (69) 

By  subiftiiuting  <66>  through  (69)  into  (30),  (41),  (49)  and  (S7),  the  function  /  occurring  in  (63) 
can  be  easily  evaluated  in  suppoit  of  the  explicit  scheme.  We  have  that 


(70) 


Coniisting  of  (65)  and  the  above  formula,  the  explicit  solver  li  easily  programmed,  yet  the  time 
step  A/  must  be  limited  to  stabilize  the  scheme  Experience  has  shown  that  we  can  simply 
enforce  the  Von  Neumann  criteria  [4J  for  the  single  step  algorithm  i.e.. 


lanilfj,) 


1 

'^2 


(71) 


In  the  context  of  this  restriction,  /,  and  are  the  lengths  of  the  cell  sides  existing  in  the  finite 

volume  grid.  It  follows  that  At  must  shrink  in  proportion  with  the  refinement  of  the  grid  in  order 
to  mainiain  the  stability  of  die  numerical  solution. 

U  Implicit  Time  Integration 

in  Its  most  fundamentai  interpretation,  impiicit  aigonthms  do  not  increase  die  accuracy  of 
a  numerical  solution  Instead,  ail  of  the  additional  macliinery  needed  by  the  implicit  scheme 
simply  permits  the  use  of  larger  time  steps  by  limiting  the  growth  of  numerical  truncation  and 
round-off  errors  By  itself,  an  implicit  scheme  can  tolerate  a  great  deal  of  error,  yet  it  does  not 
allow  this  error  to  grow.  This  concept  is  well  illustrated  by  the  work  on  model  problems 
performed  by  Risher.|l  I]  On  the  computer,  implicit  algorithms  are  more  difficult  to  program 
because  they  usually  involve  solving  systems  of  linear  equations  Linear  systems  characterized 
by  matrices  possessing  more  than  three  diagonals  are  very  difficult  and  computationally 
expensive  to  solve.  Con>;ider  that  a  typical  implicit  solver  for  the  ID  heat  equation  usually 
involves  a  tridiagonal  matnx.  For  each  additional  dimension,  we  can  count  on  adding  at  least  two 
diagonals  to  the  matrix.  The  first  step  in  our  implicit  discretization  procedure  is  to  ensure  that  we 
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vU  only  iwo  diagonal/  per  &pace  dimension.  To  accomplish  this  goal,  we  refrain  from  using 
i^iencil  averages  for  derivaiivei;  token  parallel  lo  the  cell  sides.  Why^  Derivaiive  averaging 
reijijirei!  the  u v  of  iwo  cell  values  exterior  to  ihe  normal  stencil  We  can  eliminate  these  averagei; 
by  i^eriing 


SjA  sS  u\  sf^.wl  s^.wl  sO 

There  is  a  dixadvaniage  associated  with  not  using  «tencil  averaging.  This  technique  does  enhance 
ifcheme  accuracy,  and  as  is  shown  in  the  next  section  of  this  report,  implicit  schemes  suffer  this 
losa  of  precision  On  the  other  hand,  the  enhanced  stability  provided  by  the  implicit  scheme  does 
Militate  the  use  of  larger  time  steps  and  produces  steady  state  solutions  more  quickly 

Even  if  stencil  averages  are  not  used,  without  some  simplification,  a  solution  of  die  heal 
equation  in  three  dimensions  requires  contending  with  a  heptadiagonal  matrix,  not  an  inviting 
prospect.  Instead,  we  employ  a  process  called  factorization  to  reduce  the  solution  of  a 
pentadiagonaJ  matrix  for  the  two-dimensional  problem  to  solving  a  sequence  of  two  tridiagonal 
matrix  problems.  The  resulting  method  is  referred  to  as  an  Alternating  Direction  Implicit  (ADI) 
scheme. 

To  begin  the  derivation  of  our  ADI  scheme,  we  first  write  the  combined  heat  flux  term 


=  (72) 

where. 

The  right  side  of  (72)  has  been  expressed  in  terms  of  the  sum  of  two  difference  operators.  No 
addiuonal  approximations  have  been  made  to  obtain  this  form,  but  it  does  perform  the  important 
step  of  separating  ihe  directional  difference  terms.  Equations  (73)  and  (74)  contain  the  formulas 
for  these  differences  where  and  are  obtained  from  (29)  and  ihe  discussions  in  Section  2.3. 

These  expressions  also  contain  the  notation  for  the  index  shift  operators  <)*'  )  and  (y‘*  >.  For 
example,  die  (r‘)  operator  functions  as 

=  =■.,,  (75) 

Wid)  the  use  of  (72),  (64)  can  be  discretized  with  the  fully  implicit  form 

S‘f  -  iT,  =  4/  (4 . + 4«  )«,7’  (76) 

In  operator  form.  (76)  is  written  as 

n- 4/ (4.. +4^  )jn, (77) 
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Noie  thdi  in  (77).  ^  u  lime  level  a^l  depends  on  other  values  of  S  at  od)er  values  on  ihe 

space  siencil  etisong  at  itte  same  lime  level  Ii  is  therefore  not  possible  to  evalu^ie  knowing 

only  the  values  i7,‘, .  As  a  result,  we  must  simulianeously  solve  for  at  all  values  of  i  and  j 

on  ihe  space  grid  since  these  values  are  interrelaied  Although  it  induces  some  error 

within  the  numencaJ  approKimjiion.  let  us  rewnie  (77)  as  follows. 

(1  -  A/  A;  -  d/  +  <y  *  1  =  O'  +  C*  +  0(iV  * )  (78) 

The  product  ^mply  implies  that  the  A;  and  operators  are  applied  in  sequence.  We 

have  also  assumed  that 

Moreover,  this  term  is  an  error  that  is  generally  acceptable  given  that  the  other 

difference  terms  are  0(^)  Using  factorization,  (78)  can  be  replaced  by 

(1-A/A;l(1-A/A,lC  +A/*  A.^^cr,  +0(df*)  (80) 

where  the  (7(A/’)  has  been  carried  forward  to  represent  the  approximation  error  incurred  in  (78) 
In  the  course  of  deriving  other  implicit  time  integration  schemes,  (80)  may  be  replaced  by  other 
formulas  that  are  ei^ually  valid  Note  that  the  current  form  is  used  for  this  speafic  derivation. 
We  continue  by  writing 

(1-A/A;l(1-AfA,]irr’=5,>A/A,i^>A/®A^i^-,-AfA,fi:‘,  (81) 

Simplifying  (81)  we  obtain 

(1  -d/  I  (1  A«  1 C’  =  (1 + Af  A^ )  ",  -  (1-  A/  A..  )6>  4, 5'  (82) 

Now  let 

(1 + A/  A^ )  O'  =  |1  -  A/  A^  I K,  (83) 

Although  not  rigorously  so.  die  if, \  are  accepted  as  existing  at  a  time  level  intermediate  to  n  and 
n*l.  In  tnith.  the  S,\  are  intermediate  values  (but  not  in  time)  resulting  at  the  end  of  step  one  of 
a  two-step  approximation.  By  substituting  (83)  into  (82),  we  obtain 

|1  -  A/  A..J  [1  -  A/  Afl  I C'  =  (^  -  A'  4  •><-  -  A  •  A,  K,  (84) 

Of 

(1-A/A..1(1-A/A,1C'  =l’-<l'A.l(i7;,-A»Afl< )  (85) 

By  equaling  the  arguments  for  the  common  operator  in  (85), 
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H  -  A/  A,  ]  ff,";’  =  S,\  -  A/  4,  u,'^  +  « A/’ )  <86) 

Equjiioni  (83)  and  (86)  koniiiiuie  ihe  Douglaa-Rachfurd  ADI  icKeaie  [41  for  citrvilmiur 
coordinaii!  kyuems.  Tliis  two-aiep  linie  intc$raiiun  kCheme  is  ejiecuiMJ  in  ihe  following  order. 

Siep  1 ;  (1  -(V  A^JD,;,  =  (1 + d/ A.  JD'  (87) 

Siep  2;  (1  -  AT  A,  lu,*;'  =  u,\  -A/  A,  ( S8) 

Nuie  dui  the  »cbeme  impIiLUly  sweeps  firtii  in  ihe  direction  (Step  I ),  and  second  in  the  // 
direLtiun  (Step  2).  The  aJierauiing  of  implicit  sweep  directions  is  chaucierisik  of  ADI  Hcheraed 
Again,  it  must  Pe  empluikimJ  that  the  overull  scheme  is  looilly  occuruie  to  0(At’  > . 

2^  Boundary  Cundliloits 

No  discussion  of  techniques  for  solving  IBVPs  is  complete  without  describing,  how 
boundary  conditions  ore  implemented.  This  statement  oi  paruculariy  uue  for  numerical  methods 
since  the  programming  of  boundary  condition  algunihms  can  become  complicated.  For  tbe 
problems  solved  later  in  this  report,  two  typ^s  of  boundary  conditions  ore  addressed'  (i)  tbe 
Dinchlet  condition  (17  is  prescribed  on  the  boundary  at  all  limes),  and  (li)  the  Neumann 
condition  (ds/dn  is  prescribed  at  every  point  on  the  boundary  at  all  times).  For  beat  conduction 
problems,  the  Dinchlet  condition  sets  the  acruol  temperature  at  each  point  on  the  boundary  On 
tbe  other  hand,  the  Neumann  condition  seb  the  normal  beat  flux  (or  temperature  gradient)  at 
each  point  on  tbe  boundary.  To  illustrate  the  means  by  which  these  condiuons  ore  enforced, 
consider  that  Figure  4  contains  the  domain  of  concern.  The  "i '  coordinate  lines  form  arcs  in  thi« 


figure.  A^Mine  that  the  uiiendani  grid  lines  itre  indexed  by  r  *  I . MAX.  In  the  radial  direction 

lie  the  ”/  'coordinate  Hoes  indexed  byj^  1.  ,,,  MAX,  Let  us  enforce  a  Dinchlet  condition  ui  the 
endpoint*  of  die  /  coordinate  Uttes.  For  tbe  expbeit  lime  integration  scheme  this  condition  is 
implemented  us  follows.  At  die  conclusion  of  each  time  step,  nuke  the  following  senings. 
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u, ,  =  u 


“,.1=“,..^' 


MAJ  ~  “•"«/  '  “i  JW.U' 


=u 


(89) 

(90) 


for  rsl,  JMAX.  For  ctK  Dougld^-Rachford  ADI  scheme,  unplemenimg  ihji  condioon  is  more 
complicdied.  Tlie  Dirichlei  cottdiuon  is  for  Step  2  and  can  be  wrinen  ai 


J-V3 


(91) 


for  t  s  1 . /MAX,  j  s  2,...,  JMAX~\.  simply  reprei^enis  die  known  terms  on  ihe  right  hand 

side  of  (88).  Specific  equations  involving  boundary  values  an/e  when  J  =  2  and  J  =  /MAX- 1  The 
restrictions  prescribed  by  (89)  and  <90>  are  still  applicable  for  (91).  Tliai  is  to  say.  the  term 
^^’1  I  1  ^  known  and  hence  moved  to  die  right  hand  side  of  (91).  The 


modified  equation  becomes 


In  the  same  manner,  an  equation  can  be  obtained  for  the  case  when  j  =  JMAX~L  i.e., 


(92) 


(93) 


There  are  only  two  diagonal  terms  in  each  of  equations  (92)  and  (93)  When  these  equaooru;  are 
combined  with  (91 )  (applied  fon*  &  2, . . .,  IMAX~  1 .  y  s  3,  ,  /MAX-2),  a  tridiagonal  linear  system 
1%  obtained  for  Step  2.  This  ADI  step  can  be  easily  solved  by  using  the  Thomas  Algorithm 
described  in  Appendix  A. 

The  Neumann  boundary  condition  is  applied  along  the  end  rays  of  the  annului;  shown  in 
Figure  4.  For  the  problems  solved  in  the  next  section,  we  designate  this  boundary  as  an  insulator. 
That  IS  to  say,  the  heat  flux  across  this  boundary  is  zero.  This  condition  is  enforced  by  setting 


u:*’  = 


Q 


tMtX  I 


rvo.ii  > 


“v,=-8., 


(94) 

(95) 


fory  •  2,...,  /MAX-1.  Step  I  of  Ihe  ADI  scheme  requires  this  boundary  condition:  its  difference 
equation  is  written  as 


for  I  B  2,...,  /MAX- 1,  /  s  2,...,  /MAX- 1  R'  ^  simply  represents  the  known  terms  existing  on  the 
right  hand  side  of  the  equation.  Consider  the  case  where  i  =  2;  the  specific  difference  equation  is 
wrinen  as 
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(97) 


To  iraniform  (97)  inio  a  pfoper  boutvjary  equaiion,  we  enforce  ctte  rightmo&i  equation  in  (94).  (n 
panicular.  (T, ,  iake«  ctte  place  of  ,  but  in  ihn  ca&e,  diere  is  no  fni^auon  of  lemvi  to  the  n^i 
hand  ^de  of  the  equation.  After  simplifying,  we  obtain 


bjV 


•tj 


(97) 


With  the  uie  of  (94),  a  aimilar  equation  can  be  derived  for  the  case  where  i  s  IMAX~\.  To  be 
brief,  we  simply  state  the  result. 


A/ ft 


'>  Ifktf  .(.V8 


“ia«-2,+(^-A/ft 


)ffn 


sR 


mux-',/ 


(98) 


The  system  formed  by  (96)  for)  -  3....,  /MAX-2  and  j  *  2,...,yMAX-l  istridiagonal  in  .structure 
and  may  be  solved  by  the  method  described  in  Appendix  A. 
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3.0  RESULTS  AND  DISCUSSION 


Although  the  algorithms  pee^ented  in  Section  2  of  this  report  are  bised  upon  archival 
methods,  this  effort's  end  product  differs  from  many  conventional  heal  transfer  methods.  This 
claim  IS  particularly  true  for  the  ADI  scheme  developed  for  generaliaed  coordinates  Aji  a  result, 
our  algorithms  rei^uire  telling,  or  in  more  appropriate  terms,  validation  In  the  following 
discussioni,  the  new  scheme  iii  tested  for  three  heal  conduction  problems  cast  in  two  dimensions. 
Each  problem  tests  a  diflereni  aspect  of  the  scheme’s  numerical  performance.  Moreover,  each 
problem  possesses  an  exact  solution  that  can  be  compared  one  on  one  with  ihe  numerical 
solution. 


3.1  Rectangular  Asymmetric  Heating  Problem 


The  first  test  problem  is  designed  lo  identify  programming  errors  for  the  finite  volume 
scheme,  lui  rectangular  Cartesian  geometry  is  extremely  easy  lo  gnd.  and  the  associated 
iransformaiion  terms  and  normal  vector  components  are  easily  calculated  on  paper  for 
comparison  with  the  code  calculacoiu.  The  associated  IBVP  is  presented  below. 


au 


0<x<lCtfl, 


0<,v<100. 


uauU.y.r) 

u(x.0. 0=25:  u(x.1 00,  r )  =  1 00 


^{O,y,0  ft  Ai(100.y,r)  . 

— - - u.  - - - u 

&  9 

u{x,y.0)s0 


/>0 


{99) 


The  IBVP  (99)  is  defined  over  a  large  domain  with  Dirichlei  boundary  conditions  enforced  at 
^  tO  and  v  =  100.  Neumann  (insulator)  boundary  conditions  are  imposed  over  the  linear 
boundaries  located  ai  xsO  and  4  =  100.  For  this  problem,  we  have  generated  a  uniform  grid 
with  Atb^vsI.  Note  that  unit  spatial  step  sizes  are  frequently  used  for  “benchmarking’’ 
numerical  solvers  Remarkably,  we  can  use  a  lime  step  size  of  I  unit  and  still  maintain  algorithm 
stability.  The  maximum  soluuon  time  is  set  ai  2500  units.  Note  that  ihi.s  problem  is  completely 
dimensionless.  Since  insulators  are  set  at  1=0,1 00 ,  it  is  evident  that  the  solution  for  {99>  varies 
only  in  the  >  direction.  That  is  to  say,  the  solution  is  the  same  at  every  value  of  x.  As  a  result, 
the  exact  solution  is  given  by  die  solution  of  the  heat  ei^uaiion  in  one  dimension  with  the  same 
Dinchlei  boundary  conditions.  The  exact  solution  for  this  problem  is  derived  in  Appendix  B.  A 
time  sequence  of  numerical  solutions  is  shown  in  Figure  5  At  each  distinct  time,  the  numerical 
solution  is  plotted  against  the  corresponding  exact  solution.  As  one  may  surmise  ftom  the 
dimensions  of  the  domain,  the  finite  volume  gnd  consists  of  mere  10,000  cells,  not  a  large  mesh 
by  today's  standards.  As  Figure  5  illustrates,  there  is  excellent  agreement  between  the  numerical 
and  exact  solutions  at  each  time  level  Also,  as  is  evidenced  by  Figure  6  ai  time  1500,  the 
numerical  solution  exhibits  the  required  uniformity  in  the  x  direction.  The  temperature  contours 
(or  contours  of  u  )  essentially  over  plot  one  another  This  type  of  performance  indicates  that  the 
Onite  volume  scheme  is  operating  properly.  There  is  no  indication  of  the  presence  of  numerical 
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Figure  5:  Pniblem  1  Soluiioiu;  Ciikuliiied  by  Using  ihe  Douglu-RadiAjrd  ADI  SchCAe  fnr 
Curvilinear  Coordioates.  Numencal  Soluiiona  are  Represented  by  Solid  Lutes.  Exact  Solutions 
are  Represented  by  Symbols. 


Figure  6:  Problem  I  Solutions  Calculaied  by  Using  die  Duuglas-RdCh/brd  ADI  Scheme  for 
CurvilinearCoordinaies.  Solutions  arc  Compared  Along  Lines  of  ConiUint.i  uTime  1500 

“heat  leaJti'’  from  the  gnJ,  so  the  programming  seems  to  contain  no  flaws.  Only  minor 
differences  appear  in  the  curves  near  times  41,  82  and  97.  Rjr  these  reasons,  we  can  conclude 
that  the  numerical  scheme  performs  quite  well  for  Problem  I 

3^  RudJul  AsymuMiric  Heulioe  of  an  Auoulur  SeclIOD 

The  previous  test  case  is  designed  to  test  the  fundamental  structure  of  the  finite  volume 
scheme.  It  ^so  validates  the  basic  time  scheme,  but  only  for  the  Cartesian  gnd  where  most 
numerical  schemes  excel.  In  this  case,  we  have  designed  a  test  problem  that  will  validate  both  the 
explicit  and  implicit  schemes  lor  a  non-Caitcsian  geumeti^.  In  this  caw,  the  test  case  is  defined 
on  on  annular  sector  as  shown  in  Figure  4.  The  a2tmuinal  extent  for  the  problem  ts  only  90’ 
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tnMcUitj  of  ihe  piciureJ  180*,  yet  ihe  boundary  ;:uaditions  are  m  forced  os  iDuiiiraied.  The  inner 
radiu&  of  ibe  aaaulus  ooe  while  die  outer  radius  i£  two.  The  IBVP  is  presenied  below. 


1<f<2:  />0 

di  4  4 

vsN{r,d./> 

un-^.f)*25;  «(2.ft/)*100  (100) 

du(r.;r/4.0  _^(r.9;r/4.0 

ee 

u{r.^.0)B0 

An  extinkiOiUion  of  the  [BV7  reveals  ihiti  iIil«  problem  is  aiimuihally  symmetric.  Still,  the  riniie 
volume  uigurithm  in  prognimmed  m  <1  generiiJ  miinner.  (i  compuiee  the  solution  throughout  the 
curved  two-dimensiunal  field.  This  problem  constitutes  on  cxcelleni  lest  of  the  idgunthm  since 
we  can  eoiily  compare  solutions  alo^  annular  radial  lines  to  reveal  any  bias  associated  with  tbe 
space  scbeme.  As  in  the  case  of  Problem  1,  we  utilize  a  100  x  ItX)  gnd  with  uniform  radial  and 
angular  Increments.  The  results  are  iJlusiraied  in  Figure  7.  a  plot  of  the  explicit  and  ADI 
numericnl  solutions  against  the  exact  solution  for  tunes  0.1,  0.2  and  I  0.  In  each  case,  the  explicit 
solution  agrees  extremely  well  with  the  exact  solution,  and  the  ADI  solution  also  performs  well. 
As  was  dcscussed  earlier,  the  ADI  scheme  does  suffer  from  a  reducuon  in  accuracy  due  to:  (1)  the 
lack  of  derivative  stencil  averaging,  and  (ii)  the  ADI  factorizaiiun  error.  It  is  nonetheles.>; 
impunant  to  realize  that  the  ADI  error  is  ooi  cumulative.  The  error  does  not  increase  in  time: 
also,  the  ADI  solution  cxinverges  to  the  steady  state  solution  with  accuracy  equaling  that  of  the 
explicit  scheme. 


Figure  7:  Problem  2  Time  Sequence  of  Radial  Solution  Plots.  The  Explicit  (EXPLk  and  ADI 
(IMPLk  Numerical  Solutions  are  Compared  Against  the  Exact  (EXCT)  Solution  at  Times  0  1,02 
and  1,0. 

The  accuracy  of  the  ADI  solution  improves  as  time  advances:  ultimately,  it  plotx  on  tup  of  the 
exact  solution  A  derivation  of  die  exact  solution  for  Problem  2  is  provided  in  Appendix  C 
Figures  S  examines  the  level  of  azimuthal  symmetry  that  is  exhibited  by  the  numerical  solutions. 
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BiLse«J  upun  ihe  boundury  cundiiioiu  uf  the  IBVF,  (he  eiuct  solitiion  is  perfectly  symmetric.  Buih 
etplicii  ood  implicit  numenctil  suluiions.  possess  a  high  level  of  symmedy  lo  Figure  S,  we  ploi 
the  esplicii  solution  iit  lime  fl.l  cast  aiong  radii  Ictcaied  at  45',  67  5*.  90*.  1 12.5*  aad  135*.  There 
IS  no  discemihle  difference  between  the  sotuiions  along  these  rays.  Figure  6  also  contains  a  plot 
of  (be  explicit  saluuoas  at  uoie  0.2  along  the  same  rays.  In  truth,  there  ore  differences  between 


Figure  S:  Problem  2  Analysis  uf  Aaimuthal  Symmetry  Occurring  in  Numerical  SoIu(joil«  Plot>i 
are  Shown  for  (he  Explicit  (EXPL)  Soluuon  ai  Time  0. 1  and  for  the  ADI  (TMPL)  Solution  at  time 
0,2  In  Each  Case.  Numerical  Solutions  are  Graphed  Along  the  Annular  Rays  Located  at  45*, 
67,5*.  90*.  112.5*  and  135*. 

the  ADI  solutions,  but  (he  differences  are  very  small  and  attributable  (o  the  errors  inherent  in  the 
approximations  used  to  create  the  ADI  scheme.  The  ensuing  errors  are  greater  in  magnitude  for 
curvibnear  gnds  than  fur  Cartesian  grids  due  to  the  effect  of  metric  scaling.  Still,  both  the 
explicit  and  implicit  solutions  exhibit  a  high  degree  uf  aeimutbal  symmetry  and  thus  perform 
very  well  in  this  test. 

3.3  Rudiully  and  .Aritnudially  AsymiiMtrtc  Heating  of  uii  Annular  Soedoo 

So  far.  we  have  successlully  tested  our  finite  volume  numerical  scheme  for  accuracy  in 
lime  and  for  spatial  accuracy  on  non-Cartesian  grids.  This  final  test  is  designed  to  determine 
wbeihcr  or  not  the  scheme  can  cope  with  asymmetric  azimuthal  boundary  condiiioii.s.  This  test 
case  is  more  difficult  than  either  of  Its  predecessor.^.  In  this  case,  the  outer  radial  boundary 
condiLon  IS  permitted  to  vary  with  ^  while  the  inner  boundary  value  is  held  fixed.  The  IBVP  is 
described  below.  All  problem  parameters  are  dimeosiQoles.s. 
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(I01> 


di 

u -utr.$.i) 


u{Xfi.f)  =  V.  «(1^.0  =  2+ 


u{r.^.O)sO 


The  exjct  vieady  »Ltie  Holuiion  for  ihis  problem  is  derived  in  Appendix  D  The  exeei,  lime 
depenJeni  solitikui  lu^i  oui  been  kompuied  beeouse  or  lU  cuaplexiiy.  and  iil»o  in  die  prior  leki 
cues,  the  dnte  scheme  perfumu  very  svelL  Instead,  this  lesi  case  concentruies  on  volidiiiing  the 
computer  programming  for  steady  state  convergence.  Based  upon  ilie  results  shown  lor  the  other 
test  cinies,  the  results  show  helow  ore  lor  the  ADI  scheme.  The  expliat  method  is  expected  to 
perform  well  in  ihU  problem,  so  those  results  ore  omined  Figure  9  contains  a  sequence  of  plots 
cost  along  the  radial  Une.«  (rays)  located  ai  ^  c^uol  to  4S*,  67.5*,  90*,  1 12  5'  and  US'.  Aiong 
each  ray,  the  numerical  ADI  solution  is  compared  with  the  e.xoci  solution 


Figure  9:  Problem  5  Comparisons  of  the  Numerical  ADI  und  Exact  Solutions  Solutions  ore 
Plotted  Along  Radial  Lines  Located  at  45*,  57.5*.  90*.  1 12.5*  and  135*. 

Along  each  ray,  the  numencoJ  solution  exhibits  very  good  agreement  with  the  exact  solution 
The  poure^ii  agreement  is  perhaps  obtained  nntresi  die  hut  outer  radial  boundary.  In  this  regiuc, 
the  computer  program  is  required  to  conform  to  die  natural  boundary  temperature  gradient  across 
only  one  finite  volume  cell  It  is  expected  ilui  this  performance  may  be  improved  by  re&oing  the 
gdd  near  the  outer  boundary.  Other  aspects  of  the  exact  solution  ore  caprured  quite  well  It  is  also 
interestiug  to  view  color  contour  plots  of  a  time  sequence  of  numerical  solutions.  These  plots  ore 
shown  m  Figure  10  for  solution  limes  0.1,  0.2,  0.3  and  0.5  The  snapshot  at  lime  0.5  is  included 
since  it  represents  the  nearly  converged  steady  state  solution.  The  heat  conduction  solution 
behsvex,  not  surprisingly,  like  a  classic  dlf^slon  problem.  Near  the  hot  boundary,  the 
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The  Qniie  volume  ADI  method  produces  solituoiu  ihai  capture  ihe  une  dependent  e/fects  of  ihe 
heal  cooduction  equation  Also,  this  nlgoriiboi  ui  stable  aoJ  produces  good  spatial  accuracy, 
especially  considering  the  approximationi)  necessary  to  irapleraeni  the  implicit  scheme  for 
curvilinear  grids. 


lemperanire  changes  very  t^iuckly  driven  by  high  temperature  gradients  The  elTect  Is  shown  in 
the  &r»i  three  plots  in  Figure  10  Notice  how  quickly  the  high  temperature  field  propagates 
through  the  annular  section.  Observe  also  tbat  the  lower  temperanire  gradient  near  the  right  side 
of  tbe  aooulu>;  cauies  a  slower  temperature  responiie  io  time.  Further  consideration  of  this  figure 
reveals  that  the  temperature  contours  change  very  LtiJe  between  limes  0.3  and  0.5.  The  reason  for 
the  ''slower”  propagaiion  of  the  temperature  field  is  caused  by  tbe  physics  of  the  heat  equation. 
The  lime  dependent  term  lor  solutions  of  this  cquauon  has  an  exponential  dependency.  That  Is  to 
say.  for  a  small  time  interval,  this  term  c'auses  large  changes  m  the  solution.  Yet,  for  large  time 
mtervals.  the  time  dependent  term  follows  an  asymptote  that  approaches  zero,  so  at  ever  larger 
times,  the  time  dependent  term  has  loss  effect.  Hence,  the  temperature  field  changes  more 
slowly.  A  parallel  effect,  particularly  significant  at  early  times,  is  driven  by  the  temperature 
gradienL  High  gradients  motivate  large  local  changes  in  the  temperature  field  As  we  discussed 
earlier,  this  effect  is  evident  at  times  0.1  and  0.2.  Based  upon  this  evidence,  we  conclude  thai  our 
ADI  solver  is  validated  for  two-dimensional  heat  conduction  problems  on  curvilmear  geometries. 


Time 

□.1 


Tin># 


j  s  .  '•*»»; 


Figure  )0:  Problem  3  Color  Contour  Plots  of  the  Temperature  Field  on  on  Annular  Section.  Plots 
of  the  ADI  Numerical  Solution  are  Shown  at  Times  0.1.  0.2.  0.3  and  0.5  All  SoluLon  Parameters 
ure  Dimensionless. 
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4.0  CONCLUSIONS 


This  repoit  de&cribed  il)e  numencdl  meittodology  behind  a  new  implicit  heal 
cotvjucuon  equation  solver  based  upon  ihe  finite  volume  method  foe  curvilinear  geomeirie^.  A 
principal  strength  of  this  method  is  ihJt  ii  pemuis  the  dccuraie  enfoeoemeni  of  boundary 
condiuon^  along  curved  boundaries  without  the  need  for  complicated  interpolation  or  interface 
tracking  scheme.«.  The  application  of  the  finite  volume  method  allows  an  effective  reduction  of 
order  for  die  ^tial  discretization  v;heme  avoiding  die  need  for  discretizing  the  Laplace  operator 
in  generalized  coordinates.  Also,  since  it  is  unnecessary  to  differentiate  basis  vectors  m  die  local 
coordinate  system,  the  scheme  remains  highly  accurate  with  far  fewer  gecmeiric  calculations. 
However,  caution  must  be  used  when  applying  this  scheme  for  problems  on  like  topologies  with 
cell  shapes  that  change  drastically  within  comparatively  small  regions  of  the  grid.  Without 
adequate  mesh  resolution  in  these  locales,  error  can  accumubie  rapidly.  A  proiorypical  computer 
program  with  both  expliai  and  implicit  lime  integration  schemes  has  been  developed  for 
curvilinear  coordinates  in  iwo  dimensions  and  validated  against  a  set  of  exact  solutions  for  an 
annular  geometry.  Our  numerical  solutions  have  proven  to  be  very  accurate  when  compared  with 
the  exact  solutions.  The  explicit  solver  is  highly  accurate,  and  the  pans  of  the  ADI  algorithm 
causing  higher  error  have  been  identified  Based  upon  its  simple  finite  volume,  semi-discrete 
form,  this  algorithm  is  readily  adaptable  for  problems  set  m  three  curvilinear  dimensions  This 
algorithm  has  many  potential  applications  for  the  analysis  of  engineering  systems 
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APPENDIX  A 
The  ThninftF^  Algorithm 


The  Thomas  Algorithm  used  lo  solve  a  system  of  linear  equation/  that  be  written  in  the 

form 


M«=ff 

where  M  is  a  irididgonal  matrix.  This  matrix  has  a  form  such  as 


1102) 


Ms 


(103) 


“/.* 

^.1 

All  entries  other  than  the  a,,  b,  and  v,  are  zero.  The  index  structure,  i.e.,  i  •  2,  ..../-I,  is  used  to 
mimic  the  grid  indices  u/ed  in  the  ADI  .scheme  discussed  in  Secooo/  2.5  and  2.6.  Here,  the 
maximum  index  /  corresponds  to  either  IMAX  or  JMAX.  The  solution  vector  u  has  the  form 


W  s[Uj 

while  the  vector  of  knowns  A  has  the  form 


<I04> 


R  =  \R^,  .,R,J  (105) 

The  Thomas  Algorithm  functions  by  applying  Gaussian  Elimination  to  this  linear  system.  An 
equivalent  result  may  also  be  odiained  by  the  use  of  LD-Factorizauon  (n  either  case,  for  a 
tridiagonal  matrix  M  that  is  {i)  sufficiently  well  conditioned  and  (il)  ftee  of  zero  entries  on  the 
diagonal,  a  solution  for  the  linear  system  may  be  obtained  m  two  steps. 

Step  I  IS  referred  to  as  forward  elimination  and  is  described  by  the  following  algorithm. 
Define  x, wy,  sO  widi  *0  Fori'  -  2. .... /AfAX-1.  compute  the  following quantiiies' 


R  -a  X 


—1 


-r 


(106) 


Clearly,  x,  depends  upon  .  so  the  quantities  must  be  calculated  together  in  a  loop  Step  2  of 
this  algorithm  is  referred  to  as  backward  subsuiuiion  and  is  described  as  follows.  At  the 
conclusion  of  forward  elimination,  it  is  immediately  evident  that 


“/.I 


(107) 


The  rest  of  the  solution  is  given  by  calculating  for  i  s  A2, ...,  2, 


w,  a  X,  -f  y,  M,., 


(108) 
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APPENDIX  B 

Solution  for  Heat  Equatkiii  wiib  AAynmirtrk  Heating  lii  One  Cartesian  DUnen^ilon 


In  many  ca&es.  ihe  heat  equation  is  very  easy  lo  solve  for  one-dimensional  Caitesian 
coordinaies  by  using  separation  of  vamblei;,  but  for  the  ai^ymmeinc  heating  problem  described 
below,  thin  method  requires  modification.  Tl^is  IBVP  i^  wniien  as  follows' 

—  =  — r,  0<x<E,  />0 

ar  ^ 

er  =  M<*.r)  (109) 

u<0,r)  =  Wgx0;  Mli,/)  =  M4xO;  neneral} 

uU,0)  =  0 

The  enforcement  of  non-zero  boundary  conditions  present/;  problem.s  for  the  .separation  of 
variables  method.  That  is  to  say,  if  w  is  zero  on  a  boundary,  then  this  condition  can  be  enforced 
by  making  the  npace  dependent  pan  of  the  solution  zero  at  the  same  location.  For  a  non-zero 
value,  ihi/i  technique  will  not  wortL  Fortunately,  we  can  work  around  this  difficulty  Equation 
rl09)  involve/;  a  linear  partial  differential  equation,  so  it  accepts  superposed  solutions.  As  a 
result,  we  can  write  u(x./)  us  the  sum  of  a  purely  .space  dependent  function  u.(x)  and  a  time 
and  space  dependent  function  u(x,r},  i.e, 

uU.t)su,(*)-»-u(x.r).  OsxSi,.  rkO  (110) 

The  function  u,(z)  is  defined  to  satisfy  the  boundary  conditions  of  (109);  hence,  w(x,r>  must 
satisfy  zero  boundary  conditions  at  XsO.  L.  Both  a.<x)  and  u(z,r)  satisfy  the  partial 
differential  equation,  so  for  w,(x) ,  we  have  that 


dr* 


sO;  u.(0)s«c:  w,(C)s«t 


The  solution  for  <  1 1 1)  is  easily  calculated  as 


<lll> 


For  u(x,r) ,  we  can  now  solve  a  different  IB  VP  given  by 


(112) 


^6  _  , 
a  “**’ 

K  *W(>,0 


0<x<i,; 


t>0 


u(0,0*«(i^0a0 


(113) 


ti(x,0)s-u,(x) 

Equation  (113)  is  now  in  the  conventional  form  for  a  heat  conduction  IB  VP;  moreover,  it  is 
solvable  by  using  separation  of  variables  with  Fourier  series.  Since  its  solution  is  performed  in 
many  classical  textbooks,  we  simply  state  the  result,  le.. 
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where 


(II4> 


More  specifically  for  (112), 


(115) 


C,=  —  (tf  t  cosfuff )  -  tfo )  (116) 

At  u  evident  from  1114),  (he  function  u  clearly  decays  exponentially,  so  superposed  solution 
approaches  the  steady  state  solution,  as 
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APPtNDIX  C 

Aiimulballv  Synuiielric  Solullnu  for  Heal  Equation  wlih  Asynioeirlc  Healing  of  aii 

Annular  Section 

Ttie  geometry  for  this  problem  is  shown  in  Section  2  6,  Figure  4.  Tlie  d&sociaied  IBVP  is 
described  as 


—  =  '  a<r <h'  />0 

d> 

w=Mr.O  ai7) 

u<a,r)  =  u, ;  «(^.r)=u, 

K<r.0)=0 

Since  Uus  problem  is  aiimuihally  symmetric,  the  angular  coordinate  does  noi  appear  in  <  177)  As 
was  discuised  in  Appendix  B,  we  represent  i/(r,r)  as  die  sum  of  iiuperposed  solutions  tf.(r)  and 
u(i'.f)  with 


«<r.r)*w,(r>  +  M<r./) 

The  function  a,(r)  is  the  solution  lo  the  steady  state  boundary  value  problem 


Careful  integration  reveals  that 


d*M,(r)  ,  1  d«.(r) 


sO 


dr*  r  dr 
er,{a)aw, ;  u,(b)a«j 


u,<r)*Cjln(r)+C, 

where 


"’ln<b>-ln(a)'  ln(b)-ln(a) 

The  ume  dependent  function  u{rj)  is  a  solution  to  the  IBVF 


dt  rdr[  dri’ 

V9u{ri) 


0<r<b; 


i>0 


U{aj)au{b.l)aO 

u{r.0)a-u,lr) 

To  solve  (122),  we  apply  the  separation  of  variables  method,  i.e .  let 


u{r.r)»/?(r)T(r) 


(118) 


(119) 


(120) 


(121) 


<I22) 


<I23> 


By  substituting  <  123)  into  die  operator  in  (122),  we  obtain  the  separable  equation 
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1  i/r 


1  d‘R 


1  ilR 


R  dr  rR  dr 


(124) 


Experience  «how&  ihai  we  can  select  a  aeparaiion  constant  with  the  value  -a  -  and  we  obtain 
two  ordinary  differenual  equauon^  ibai  can  be  solved  based  upon  (he  value  of  X  •  Hence. 


d! 


d‘R  dR 


(125) 


<I26> 


Equauon  (126),  along  with  (be  auendani  boundary  condiuoni,  consutuiei;  a  Sturm-Liouville 
problem,  so  (he  only  admissible  values  of  are  real  and  positive.  1 12]  Of  course,  (126)  is 
Bessel's  zeroth  order  diflereniial  equation  with  parameter  X.  Tbe  solution  (o  (125)  is  very 
simple  (o  derive,  so  we  slate  (be  result. 

r(f)  =  -4exp(-^*/)  (127) 

where  4  is  a  constant  (hat  must  be  obtained  from  applying  (he  initial  condiuons.  Tbe  solution  of 
<I26)  IS  written  as 

«(r)=s5^1^)+Cro(^)  (128) 

where  7g  and  are  the  zeroih  order  Bessel  and  Neumann  functions.  B  and  C  are,  at  (bis  stage 
of  (be  solution,  arbitrary  constants  pending  (he  application  of  boundary  and  initial  conditions. 

Based  upon  (he  strucrure  of  (122),  (he  boundary  conditions  for  (his  problem  are  reduced 
to  solely  space  dependent  forms,  1  e.. 


Ria)sR(h)sO 

Therefore,  by  applying  (126).  we  obtain  the  rwo  simultaneous  equations 


(129) 


B7o(.Ui)+CKo(A6)aO  <I30) 

BJo(.Ah)4.CYf,(U)^0  (131) 

Tbe  only  variable  parameters  in  (130)  and  (131)  are  5,  C  and  X  .  Clearly,  B  and  C  cannot  be 
zero,  or  (he  thvial  solution  would  result.  Instead,  let  us  rewrite  (130)  and  (131)  in  matrix  form 
where 


Since  this  system  is  homogeneous  and  linear,  die  only  way  (hat  (132)  can  be  satisfied  is  if 


(132) 


/,(/U/)  Y^^^4I) 
iM)  Y,(Xh) 


(133) 
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The  ruuii  of  (133)  ltojuIsi  of  a  couiuuM)'  inflfuie  sei  coctHiuing  of  dj»iinLi  valueit  of  A.  Figure 
conieiiiK  i  ploi  of  the  deiernkifiiini  .  You  CiUi  we  wveriil  rooui  of  IhiA  ei^uaiion  udioiied  in 
the  plot. 


*1  k  '•!*  •  >  * 


Figure  1 1:  BuunJor)'  CondiiionDeiermiiuni  VuJuej  for  ihe  Asymmetric  healing  Problem  Set  oa 
un  Annular  Domain  with  Azimuthal  Symmetry, 

Each  V  al  ue  of  saii>i(y  ing  (133)  defines  an  eigenfunction  of  ( 1 28),  A  more  specific  form  of  the 
eigenfunctiOD  (128)  may  be  developed  as  follows.  For  the  root  /(,,  of  ( 1331.  a  renrrangejneoi 
of(13 1 )  shows  that 


Cs-B 


Joiib) 

YAih) 


(134) 


The  sobsiiiuiion  of  ibis  relaiioo  into  ( 12S>  renders  ihe  eigenfunction  A(^r).  i.e. 


-'eU.o  nu.'-) 


j^a.b)  Y^a.b) 

Since  B  is  JUKI  un  arbitrary  consuinL  we  can  write  this  eigenfunclion  as 


(135) 


:i36) 


/oU,6)  Y^iXb) 

Eigenfunctions  (136)  ore  now  multiplied  by  the  solution  for  TU)  and  expanded  os  Fourier  series 
solving  u(r,r).  By  using  (118).  ;12Q)  and  (1211.  the  general  solution  for  (117)  is 


yo(V> 


uM) 


(137) 


subject  to  the  loiual  condition  from  (117) 


ii(r,0)s0 


(138) 
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Recall  (hai  Bes&el  ftinctions  (and  hence,  linear  con^buutions,  e  § .  <  I36>  of  Bessel  functions)  are 
oitttogonal  under  ihe  inner  product  [12] 

* 

Jr /Ha)  Jr  <139) 

# 

Therefore,  when  we  evaluate  <I37>  at  the  uuiial  condition  (13S)  and  form  die  inner  integral 
{u<r,0X  R(A^)).  we  can  solve  for  the  coefricienis  B. ,  i.e., 


.  _  Jr«,(r)f?(4)Jr 

JrR*(VW^ 


<I40> 


In  practice,  integralji  involving  Bessel  functioru;  can  be  time  consuming  lo  evaluate,  so  for  the 
results  shown  in  this  report,  the  integrals  in  (140)  are  numerical!)'  evaluated  b)  using  the 
trapezoidal  nile  over  some  10,000  subintervals  |7]  This  process  has  been  demonstrated  to 
produce  very  accurate  solutions. 
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APPtNDIX  D 

FuJIy  Asvnmieirlc  Heating  Problem  for  an  Annular  Section 


This  problem  jiolved  ia  order  to  provide  a  le&i  for  our  finite  volume  heat  conduction 
code  that  incorporaiet  treating  asymmetry  m  two  separate  djreciions.  Thit  problem  is  defined  on 
a  cylindrical  geomei/y,  so  the  unsteady  solution  involves  complicated  Bessel  functions.  Since  the 
temporal  performance  of  the  computer  program  has  been  validated  on  two  other  problems,  the 
exact  solution  that  follows  addre.«ses  a  steady  state  temperature  profile  an  annular  section.  The 
boundary  conditions  incorporate  the  desired  asymmetry.  Tlie  BVP  is  defined  as  follows. 


7*«a0.  fl<r<i>. 
usuir.ff} 

du{r,$^)  _du{r.&g) 

dd  ’  d0  ' 

The  Laplace  equation  above  can  be  expanded  in  cylindrical  coordinates  as 

1  *  1  . 
dr  r  dr  f  dB 


(141) 


(142) 


If  we  assume  that  .  ( 142)  is  separable  producing  the  two  ordinary  differential 

equations,  i.e.. 


dr  dr 


(143) 


£!2+^*0aO  <I44) 

dr 

for  a  separation  constant  ^ .  The  process  of  solving  (143)  and  (144)  is  complicated  by  the  fact 
that  zero  is  an  admi.ssible  value  for  .  As  a  result,  each  ordinary  differential  equation  has  two 
families  of  solutions,  one  for  ^  sO,  die  other  for  ^  *  0 . 


We  begin  by  considering  the  radial  equation  (143).  Suppose  that  ^sO:  then  (143) 
becomes 


This  equation  has  die  solution 


d^R  rfff  ^ 


<I45> 


^(r)*  Aln(r)*B  <I46) 

where  A  and  B  are  arbitrary  constants  selected  based  on  the  boundary  conditions  Suppose  now 
that  2  0 .  Equation  (143)  can  be  solved  by  assuming  a  solution  of  the  form  .  where 

a  IS  a  parameter  based  upon  Substitution  of  r"  into  (143)  yields  the  expression 
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(147) 


For  d  non-trividi  &oluiion,  is  non-zero,  so 

^(a-D  +  G-A^aO  ^  O-atA  (148) 

TtKrefore,  iwo  possible  solutions  result  for  A  s  0 ,  so  iheir  combinauon  forms  n  solution  for 
<l43).i,e.. 


R^{r)  =  Ar*+Br  *  <I49) 

Now  consider  the  angular  equation  (144).  The  case  where  AbO  constitutes  a  perfecil)' 
acceptable  solution,  i.e . 


©<,(^)aCtf  +  I> 


<I50) 


where  C  and  D  are  arbitrary  constants.  On  the  other  hand,  suppose  that  A  sO  The  solution  for 
<  144)  IS  well  known  for  this  circumstance.  In  fact. 


Qf^sCcosU^+OsiiXA^  (151) 

where  C  and  D  are  again,  arbitrary  constants.  Using  die  angular  solution  forms  given  In  (ISO) 
and  (151),  we  may  begin  to  evaluate  some  of  the  arbitrary  constants.  This  process  must  be 
accomplished  with  care  because  arbitrary  constants  C  and  D  have  different  values  for  the 
A  s  0  and  Z-tO  solutions.  Recall  die  angular  boundary  condition  ei^uations  from  ( 141  > 


dS  d^ 


Similarly, 


d$ 


d0  d$ 


dQ{$t) 


d$ 


=0 


By  applying  ( 1 52)  to  ( 1 50),  we  obtai  n 

^  =  C=0  ^  c=o 

d$ 

The  application  of  (153)  confirms  this  results:  therefore 

Forthecaiie  where  A  sO,  let  us  apply  (152)  to  (151).  We  obtain 

de(^) 

Hence. 


d$ 


= -A  CsiiX  A  )+ A  Z>cosa  )  =  0 


(152) 

(153) 

(154) 

(155) 

(156) 
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(1S7) 


cosUft) 

Equujon  (151)  can  now  be  rewniien  as 


e{e)’>C[cosJAe)co6tAei)+iaiXe>iiriXei)]~CccsAi^-$,}  (158) 

In  (158),  (he  consiani  terms  appearing  in  {157)  have  ^mply  bean  rolled  up  mio  a  single  arbitrary 
constant  C.  Now  we  apply  (153)  to  (158).  and 

C.^sO  =>  «in^(^-64)=0  (159) 

Equation  (159)  cleariy  allows  the  non-26ro  eij^nvalues  to  be  idencifial  as 


(160) 

Finally,  the  angular  eigenfunction  for  ^  s  0  can  be  wnrien  as 

0.(^)sCcos^(^-^^).  ns1,2..  (161) 

Tbe  next  step  in  the  solution  procedure  involves  building  a  general  solution  that  satisfies 
ibe  radial  boundary  conditions  given  in  (141).  Note  (hat  In  general,  (be  boundary  values  at  rsa 
and  'sh  ift  non-tnvial  In  order  to  enforce  non-tnvial  boundary  values  on  both  of  tliese  radial 
boundaries,  we  must  use  linearity.  That  is  (o  say,  u(r,r)  must  be  composed  of  two  solutions 
u,ir.e)  plus  such  that 


u{r.^sW,(r.i9)  +  Wj{r.^ 

K,(o.^s«,(^):  ti,(^,i9)s0  (162) 

If  the  components  solutions  u,  and  u,  ^  suniciured  in  accordance  with  (162),  (he  general 
solution  u(r,0)  satisfies  tbe  required  boundary  conditions.  Let  us  determine  u,.  For  ^sO,  we 
apply  (146)  and  (155)  to  obtain 

u,lh,0^siAMh)+B)DsO:  D*0  9  £s-/4Mb)  (163) 

Hence,  for  b  0,  we  have  (hat 

H,  (r.i9,  ^  s  0)  s  4ft  (Ini  r)  -  ln(b))  { 1 64} 

where  represents  a  roll-up  of  arbitrary  constants.  The  remaining  part  of  u,  is  provided  by  an 
analysis  of  die  case  where  A  sO  Apply  (149),  (158)  and  (162),  and  we  can  write 

u,(r.e.A^Q)  =  iAb^+Bh'^).CooiUe-e^)=0  (165) 

In  general,  the  angular  pan  of  (165)  is  non-zero,  so  we  conclude  that 
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<I66> 


Ab^  +  Bb-'^sO  ^  Bs-Ab^ 
Ttterefore,  with  ihe  of  iinoilier  rolled  up  con&idnt  A, , 


u,{r.eA^O)  =  A.{r'b'''-r  <I67) 

Wid)  the  u&e  of  (160),  (164)  aiid  (167),  we  c^in  write  the  general  jiolution  for  u,,  ie.. 


«,(r,^)s4o(ln<r)-JiKh))+j;A.(rn-*'-r-«fe*')cos.l.(^-^)  (168) 

A&1 

Evpan&ion  coefTicienis  A,  must  be  determined  based  upon  the  boundary  condiiioni;  enforced  at 
r-u .  For  ihii;  purpose,  we  exploit  ihe  orthogonality  of  eigenfunctiom;.  In  particular,  it  can  be 
i;hown  that  the  set 


2«(1,C08^(^-^^),n»1,2,...l 

IS  orthogonal  under  the  inner  product  defined  by 


(169) 


(170) 


By  computing  inner  products  carefully,  it  can  be  shown  that 


1 


v4.s 


A^a - 1 - f'u  W)dfi 


df) 


(171) 

(172) 


Equations  (170)  through  (172)  form  a  general  ^solution  for  u,  that  satisfies  the  radial  boundary 
conditions  u,{a.6)»  u^{&)  and  u,ib,&)s0. 


To  complete  die  general  jiolution  for  (141).  we  also  must  construct  the  component 
solution  that  conforms  to  the  boundary  conditional  given  in  (162)  Tlie  same  procedure 

used  to  obtain  w,  i.<  applied  in  this  case.  Since  the  procedure  is  repetitive,  only  key  results 
encountered  in  the  course  of  the  derivation  are  presented.  By  noting  that 
components  of  u;  may  be  obtained  as 

As  0)  S  A»(ln(r)-  lixn))  (173) 

u,(r,i?.A*0)sj4.(i/**r*-o*r-*')cOS4-(^-^)  (174) 

Ordiogonality,  (I69)  and  (170),  can  be  used  in  conjunction  with  the  boundary  condition 
iif{b,ff)su^{ff)  to  determine  coefficients  A, .  With  some  careful  calculation,  we  obtain 
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(J75) 


computer  can  be  lo  calculate  u{r,$)  to  any  desired  level  of  dccurdcy. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS  AND  ACRONYMS 


ADI 

BVP 

CEM 

CFD 

IBVP 


Aliemaiing  Direction  Implicit 
Bound^iry  V^ilue  Problem 
Compuidiional  Elecironugneiics 
CompiiiaiionaJ  Ruid  Dynanucs 
Imtiiil  -  Boonddiy  Value  Problem 
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